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Abstract 



The effect of the Berendsen thermostat on the dynamical properties of bulk 
SPC/E water is tested by generating power spectra associated with fluc- 
tuations in various observables. The Berendsen thermostat is found to be 
very effective in preserving temporal correlations in fluctuations of tagged 
particle quantities over a very wide range of frequencies. Even correlations 
in fluctuations of global properties, such as the total potential energy, are 
well-preserved for time periods shorter than the thermostat time constant. 
Deviations in dynamical behaviour from the microcanonical limit do not, 
however, always decrease smoothly with increasing values of the thermostat 
time constant but may be somewhat larger for some intermediate values of 
tb, specially in the supercooled regime, which are similar to time scales for 
slow relaxation processes in bulk water. 
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The ideal ensemble to extract dynamical information from a molecular 
dynamics simulation is the microcanonical (NVE) ensemble ^HJ. Since the 
total energy is conserved in this ensemble, the Newtonian equations of mo- 
tion can be assumed to represent the natural evolution of the system, subject 
to the accuracy of using classical mechanics to describe the atomic dynam- 
ics. The constant energy (E) and volume (V) conditions do not, however, 
correspond to the most common experimental conditions and therefore it is 
often desirable to implement MD simulations in the canonical (NVT) and 
isothermal-isobaric (NPT) ensembles. An additional reason for preferring 
the NVT ensemble is that when long run lengths are required, as in the case 
of studies of slow dynamics in liquids or glasses, there may be significant 
energy drift in the NVE ensemble. 

Two types of approaches have been developed to adapt MD simulations 
to the canonical ensemble: (i) extended Lagrangian methods, such as Nose- 
Hoover thermostats and (ii) resampling or rescaling of velocities, as in the 
case of the Andersen and Berendsen thermostats. The extended Lagrangian 
methods will generate the true canonical distribution of velocities. The 
rescaling approaches will only ensure that the average kinetic energy of the 
system corresponds to the expected value at the desired temperature but 
have the advantage that they can be combined very simply with the Verlet 
algorithm. In both approaches, the degree of perturbation of the real time 
evolution of the system can be adjusted by manipulating various thermostat 
parameters. 

The Berendsen thermostat represents a proportional scaling of the ve- 
locities per time step in the algorithm with the scaling factor being given 

3 



by 



At / To 



) 



(1) 



where At is the time step, tb is the time constant of the Berendsen thermo- 
stat , To is the desired temperature and T is the instantaneous temperature 
PJ.The velocity rescaling can be incorporated easily into the Verlet leapfrog 
algorithm. By varying the thermostat time constant, tb, one can, in ef- 
fect, increase or decrease the degree of coupling to an external bath. The 
limit, Tb oo, represents the microcanonical ensemble. The original paper 
of Berendsen et al applied this thermostatting procedure to simulations of 
216 water molecules bound by the rigid SPC potential. Static structural 
averages were shown to be unaffected for values of tb as small as 0.01 ps. 
Fluctutations in global properties were strongly affected when tb was less 
than O.lps implying that analysis of fluctuations cannot be used to deter- 
mine observable properties. Single-particle properties, including dynamical 
quantities such as the diffusivity and orientational correlation times, were 
found to agree, within the limits of statistical error, with the results derived 
from the microcanonical ensemble for tb > 0.1 ps. While the Berendsen 
thermostat is widely used, we have come across only a limited set of sub- 
sequent studies in the literature which inspect in detail the nature of the 
associated trajectories jUISlini- 

In this note, we perform some detailed tests of the effect of the Berendsen 
thermostat on the dynamical properties of water. The motivation for these 
tests comes from prior work on molecular dynamics simulations of bulk water 
where power spectra associated with a number of dynamical quantitites were 
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shown to possess a 1//° frequency regime, indicative of multiple time scale 
behaviour due to hydrogen bond network re-organisations [71 |H1 1^1 1101 lllj . 
The power spectrum is defined as 



where A{t) is any time-dependent mechanical quantity and (A) is the corre- 
sponding average over the system trajectory. Most recently we have exam- 
ined power spectra generated from fluctuations in tagged particle quantities 
sensitive to the local molecular environment, such as the tagged particle 
potential and kinetic energies, over a wide temperature range Fluctua- 
tions in the tagged particle potential energies, but not in the kinetic energies, 
were shown to give rise to 1 /f" noise, indicative of a structural origin for the 
multiple time-scale behaviour; moreover, the exponents and frequency range 
of 1//" behaviour were shown to be temperature-dependent. This suggests 
that the 1//" behaviour of power spectra may provide a simple and conve- 
nient way of monitoring changes in the hydrogen bond network dynamics 
with temperature or density. Since many of the anomalous features of wa- 
ter, originating from its unusual network structure, are more marked in the 
supercooled regime, it is of interest to study the 1//° behaviour at low tem- 
peratures. Molecular dynamics simulations of supercooled water typically 
require long run lengths of the order of 2ns or more if reliable estimates 
of dynamical quantities, such as the diffusivity, are to be obtained. The 
Berendsen thermostat is frequently used to generate the trajectory for such 
long runs to avoid the energy drift associated with micro canonical runs and 
facilitate comparison with experimental data |12| ll^^ l ITU I15 | To check 





5 



the effect of the Berendsen thermostat on temporal correlations associated 
with fluctuations in tagged particle quantities, we have carried out a set of 
tests on a range of state points for SPC/E water. 

The potential energy surface of water was represented by the rigid, non- 
polarizable, three-site SPC/E model PJj. The molecular dynamics (MD) 
simulations were performed using the DL_POLY software package pi llHirTO] . 
A cubic simulation cell containing 256 SPC/E water molecules was used. 
Electrostatic interactions were evaluated using the Ewald sum approach. 
The MD trajectory was propagated in the micro canonical (NVE) ensemble 
using the Verlet leapfrog algorithm in conjunction with the SHAKE algo- 
rithm to implement bond constraints HI. MD trajectories were also gener- 
ated in the NVT ensemble using the Berendsen thermostat with a range of 
time constants. In keeping with earlier work, a time step of 1 fs was used 
for all the simulations and production run lengths were kept at 2ns. The 
coupling constant for the Berendsen thermostat was assigned values of 1, 
10, 25, 50 and 200ps. Table 1 summarises our results from from the NVT 
and NVE simulations at the following temperatures along the 1 g cm"'^ iso- 
chore: 229K, 260K and 295K. The melting point of SPC/E water at 1 atm 
pressure is known to be below 260 K PU]. The pressure, P, and the config- 
urational energy, U, correspond to simple thermodynamic averages and the 
associated errors were estimated by block averaging 21 . The large error 
bars on the pressure are typical of systems with electrostatic interactions; 
within statistical error, our results agree well with those in the literature. 

We now consider the power spectra generated from fluctuations in the 
tagged particle potential energies at the different state points (see Figure 
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1). The tagged particle potential energy, u{t), is defined as the interac- 
tion energy of an individual molecule with all the other molecules in the 
system and the corresponding power spectrum is denoted by Su{f)- Since 
the potential energy surface is assumed to be pair- additive, the total poten- 
tial energy, U{t) = 0.5 J2iUi{t) where the sum extends over all molecules. 
Tagged particle potential energies of 32 molecules in each simulation were 
sampled at intervals of lOfs and Fourier transformed to generate the power 
spectra using standard FFT routines [12 • When generating power spectra, 
a square time window of 328 ps width and a sampling interval of 10 fs was 
used. This choice of sampling interval corresponds to a Nyquist frequency 
of 1666 cm^^ and sets an upper bound on the frequency range over which 
we can obtain power spectra. The values of tb provide the lower limit on 
the frequency range over which we can obtain reliable power spectra; thus, 
Tb = Ips corresponds to a lower frequency limit of 33 cm~^. The normali- 
sation convention was chosen such that the integrated area under the S{f) 
curve equalled the mean square amplitude of the time signal. Statistical 
noise in the power spectra was reduced by averaging over overlapping time 
signal windows as well as over individual tagged particle spectra; typically 
11 windows were used to cover the 2ns run length. In a given frequency 
interval showing 1//° behaviour, linear least squares fitting of lnS'(/) was 
done to obtain the a values with an estimated standard error of less than 
3%. Additional computational details regarding calculation of power spectra 
are given in ref.|llj. 

Figure 1(a) shows the Suif) spectra at 230K obtained from the NVE 
run as well as from NVT runs with different values of tb are shown. At this 
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temperature, the Su{f) curve shows three distinct features: (i) a broad peak, 
originating from Hbrational modes, centred at approximately 600 cm^^; (ii) 
a 1//" regime between 60 and 298 cm~^ and (iii) a second, low-frequency 
1//"^' region between 1 and 40 cm~^ such that a a'. On a logarithmic 
scale, the power spectra obtained from the different simulations are indis- 
tinguishable. Careful inspection of the >S'„(/) curves, however, shows that 
the noise in the results obtained with tb = 25ps is somewhat larger than 
in all the other cases. To quantify the effect on the Su{f) curves, we have 
tabulated the a and a' values in Table I. The low frequency exponent a' 
from all the runs agrees to within 3% except at tb = 25ps. The spread in 
the values of the high-frequency exponent a is somewhat larger with again 
a maximum deviation at tb = 25ps. While the deviation at tb = 25ps is 
not large, it is interesting that it should be greater than the deviation seen 
for the smallest tb value of Ips. The results for 260K, shown in Figure 
1(b) and Table I, are very similar in that variations with tb values are very 
small and within the estimated statistical error. The a and a' values are, 
however, very similar, indicating that it may be better to characterise the 
region from 1 to 300 cm~^ as a single 1//°' region. The curve at 295K 

clearly shows only one 1//"^ region which merges into the librational band 
with a crossover to white noise behaviour below 2 cm~^. 

In addition to the results shown for the tagged particle potential en- 
ergy, we have checked the convergence, with respect to tb, of power spectra 
associated with a number of other tagged particle quantities, such as the 
local kinetic energies and order parameters, and found similar trends to 
those discussed above. From our results, it is evident that the Berendsen 
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thermostat preserves the temporal correlations in tagged particle quantities 
remarkably well over a wide frequency range, even for Berendsen thermostat 
coupling constants as small as 1 ps but the small discrepancy seen at 230K 
for tb = 25ps does remain. Figure 2 shows the mean square displacement 
as a function of time at this state point, as an additional check, from the 
different simulations. On a logarithmic scale, the curves are indistinguish- 
able but it can be seen that the diffusivity estimates obtained at tb = 25ps 
would be slightly different than for all the other cases. 

To understand the discrepancy seen for tb = 25ps at 230K, we have 
performed simulations for three additional values of tb at 230K and 1 g 
cm^^; all the tagged particle power spectra are shown in Figure 3. It can 
be seen that the tb values of 25ps and 20ps result in significantly greater 
noise and change in shape of the Su{f) curves in the 100 to 400 cm~^ 
regime. These values of tb must therefore correspond to some physically 
significant dynamical time scales of bulk water. The frequencies of the 0-0 
stretching and 0-0-0 bending modes are at 200 cm^^ (0.165 ps) and 50 
cm~^ (0.66ps) respectively (23] • Clearly all the tb values chosen by us lie 
above the time scales associated with these intermolecular vibrations. In 
reviewing the literature on the dynamics of supercooled SPC/E water, we 
found, however, that the time at which the non-Gaussianity parameter is a 
maximum is 20.6 ps at 224K and 15.5ps at 238K |23]. The time at which the 
non-Gaussianity parameter peaks is an indication of the time scale on which 
dynamical correlations are most pronounced. This suggests that anomalous 
noise effect seen for tb values of 20 and 25ps is due to a coincidence of 
the time scales for the thermostat and for correlated motion in supercooled 
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water. At higher temperatures, such long-hved correlations are much less 
pronounced and therefore variations in the power spectra with tb are not 
as noticeable. 

Figure 4 shows the power spectra obtained from fluctuations in the to- 
tal potential energy of the system from NVE and NVT simulations. As 
expected, these power spectra are much more noisy than those obtained 
from tagged particle potential energies. Despite the fact that fluctuations 
in global quantities are expected to be much more subject to distortion on 
thermostatting than single-particle quantities, the power spectra obtained 
from the NVT ensemble simulations with tb = 200ps agree well with the 
NVE results for a frequency range of 0.2 to 1000 cm~-^. The results for 
Tb = Ips, do, however, show large deviations in the low frequency regime 
below 10 cm~^. 

This study indicates that NVT ensemble simulations using the Berend- 
sen thermostat are successfully able to reproduce the temporal correlations 
in tagged particle quantities over a very wide frequency range provided suf- 
ficiently large values of the Berendsen time constant are used. Even the 
power spectra associated with global quantities, such as the total potential 
energy, are well reproduced for frequencies greater than 1/tb- It is, how- 
ever, worthwhile to check the dynamical behaviour for a range of thermostat 
time constants since some values may show slightly larger deviations from 
microcanonical behaviour than others. This effect is more pronounced on 
supercooling and is noticeable for values of the Berendsen thermostat time 
constant which are similar to the time at which the non-Gaussianity pa- 
rameter is maximum, suggesting that the anomalous values of tb fall in the 
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range of time scales associated with slow structural relaxations. Though the 
Berendsen thermostat is widely used, we are not aware of any previous study 
which has demonstrated this effect. We expect that analogous deviations 
from microcanonical behaviour will be seen with other thermostats when 
slow relaxation processes are present. 
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Figure Captions 



1. Comparison of power spectra of the tagged particle potential energy 
fluctuations from the NVE ensemble and from NVT runs using dif- 
ferent values of the Berendsen coupling constant along the 1 g cm~^ 
isochore at (a) 230K, (b) 260K and (c) 295K. Insets show comparison 
over a relatively small frequency range between the NVE results and 
the results for a single value of tb- The figure key is the same in all 
cases. 

2. Mean square displacement (MSD), in units of A°^, as a function of 
time (in ps) for different values of the Berendsen coupling constant at 
230K and 1 g cm~^. Inset shows the data using the logarithmic scale. 

3. Comparison of power spectra of the tagged particle potential energy 
fluctuations from the NVE ensemble and from NVT runs using difl'er- 
ent values of the Berendsen coupling constant at 1 g cm~^ and 230K. 
The power spectra for different values of tb are multiplied by different 
factors to facilitate comparison of the overall shapes of the curves. 

4. Comparison of the power spectra of the total potential energy fluc- 
tuations at 295K and 1 g/cm^ for different values of the Berendsen 
coupling constant. 
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Table 1: Summary of simulation results from microcanonical(NVE) and 
canonical (NVT) ensemble simulations of SPC/E water at state points along 
the 1 g cm~^ isochore. Each run is 2ns long. The error bars on the temper- 
ature (T) apply in the case of the micro canonical (NVE) ensemble run. The 
exponent a is obtained by fitting £'«(/) to a 1//" form over the frequency 
ranges 60-298 cm-^ (230K), 60-298 cm'^ (260K) and 4-200 (295K). The 
exponent a' is obtained by fitting to a 1//"^ form in the frequency range 
l-40cm-i at 230K and 260K. 
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1.31 


1.16 




25 


-23±4 


-51.2±0.4* 
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Fig 3 
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